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azimuthal-shear  waves  (SH).  Likewise,  incident  SV  and  SH  waves  both  generate  radial  P  waves,  and  SV 
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stants  of  the  anisotropic  inclusion.  These  results  may  help  explain  on  a  theoretical  basis  the  observed 
short-period  SII  and  SV  waves  from  underground  explosions  at  teleseismic  distances. 
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ABSTRACT 


In  the  absence  of  local  scattering  the  ground  motion  due  to  major  phases  at  an  array  should 
be  perfectly  coherent  between  different  seismometers.  We  have  studied  the  coherency  of 
ground  motion  for  regional  phases  such  as  Lg  as  a  function  of  frequency  and  spatial  separat  ion 
for  the  NORESS,  ARCESS  and  FINESA  arrays.  Events  examined  are  quarry  blasts  at  a 
distance  of  200-400  km  from  the  array.  Coherency  was  estimated  for  10-25  s  windows 
containing  the  phase  and  with  the  array  seismograms  time  shifted  to  remove  the  average  effect 
of  propagation  across  the  array.  In  the  1-10  Hz  range  coherency  decreases  with  increasing 
spatial  separation.  The  decrease  is  faster  for  higher  frequencies,  but  if  the  separation  is  scaled 
to  the  wavelength  then  the  decay  curves  are  similar  and  indicate  that  coherency  decreases 
to  less  than  0.5  within  about  a  wavelength.  To  study  this  problem  we  have  simulated 
synthetic  seismograms  by  two  methods.  Finite  difference  calculations  have  been  made  for 
a  source  below  a  layer  over  a  half  space  to  produce  seismograms  for  an  array  at  the  top  of 
the  layer.  The  layer  and  the  half  space  contain  random  velocity  variations  and  the  layer — 
half  space  interface  is  randomly  rough.  Models  using  self-similar  autocorrelations  for  the 
velocity  and  interface  perturbations  reproduce  the  major  features  of  the  observed  coherency 
relations.  Synthetic  seismograms  have  also  been  constructed  by  taking  an  observed  trace 
and  producing  new  traces  by  successively  randomly  perturbing  the  Fourier  amplitude  and/or 
phase.  Perturbations  of  phase  (or  equivalently  travel  time)  have  a  much  greater  effect  on 
coherency  than  perturbations  of  the  amplitude. 
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INTRODUCTION 


At  regional  distances  (100-2000  km)  seismic  energy  from  shallow  events  travels  through  the 
lithosphere.  While  traditional  layered  models,  including  large  scale  lateral  heterogeneity, 
have  succeeded  in  explaining  travel  times  and  the  gross  amplitudes  of  regional  phases,  it  is 
clear  that  scattering  has  an  important  effect.  Three  effects  are  notable.  One  is  the  presence 
of  coda,  scattered  energy  arriving  at  times  that  cannot  be  explained  by  layered  structures 
(Aki  and  Chouet,  1975;  Herrmann,  1980).  A  second  is  attenuation  due  to  scattering  (Dainty, 
1981,  1984;  Wu,  1982;  Toksoz  et  a/.,  1989).  Both  of  these  effects  are  primarily  controlled  by 
backseat tering.  The  third  effect  is  the  incoherence  of  regional  seismic  phases  even  at  closely 
spaced  stations,  often  called  fluctuations.  This  effect  is  controlled  by  forward  scattering 
(Chernov,  1960).  In  this  paper  we  will  examine  this  effect  on  regional  seismograms.  While 
some  work  has  been  done  on  regional  phases  such  as  Lg  (Der  ct  al. ,  1984;  Ingate  et  al., 
1985),  the  effect  has  been  studied  extensively  for  teleseismic  arrivals  (Aki,  1973;  Capon. 
1974;  Flatte  and  Wu,  1988),  and  some  work  has  also  been  carried  out  in  the  field  of  strong 
motion  studies  (Farjoodi  et  al.,  1985;  Loh,  1985;  Harichandran  and  Vanmarcke,  1986;  Toksoz 
et  al. ,  1989).  Our  goal  is  to  measure  the  coherence  of  the  regional  phase  Lg  in  Fennoscandia 
and  offer  a  preliminary  interpretation  in  terms  of  the  features  of  the  wavefield  and  the  earth 
that  might  be  causing  it.  All  of  the  areas  studied  have  Precambrian  or  Paleozoic  crystalline 
rock  exposed  at  the  surface. 

Regional  arrays  such  as  NORESS  (Bungum  et  al.,  1985),  FINESA  (Korhonen  et  al.. 
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1987)  and  ARCESS  (Mykkeltveit  et  al. ,  1987)  provide  powerful  tools  to  examine  the  effect 
of  scattering  on  seismograms.  In  previous  investigations  at  NORESS  (Dainty,  1985,  1989; 
Dainty  and  Harris,  1989)  we  have  looked  at  coda,  which  consists  of  scattered  energy  often 
deviating  at  large  angles  from  the  direct  arrivals.  In  such  a  situation  frequency- wavenumber 
analysis  has  proven  to  be  useful  because  the  arrays  can  resolve  the  differences  in  angle 
clearly.  When  the  effect  of  scattering  on  direct  phases  such  as  Lg  is  considered,  however, 
this  method  does  not  produce  useful  results  because  the  angle  of  deviation  of  the  (forward) 
scattered  energy  is  not  resolvable  by  the  array. 

Previous  work  on  this  problem  has  used  the  fluctuations  of  amplitude  and  phase  (arrival 
time)  of  teleseismic  P  at  the  LASA  (Aki,  1973;  Capon,  1974)  and  NORSAR  (Capon  and 
Berteussen,  1974;  Flatte  and  Wu,  1988)  arrays  to  estimate  statistical  parameters  of  P  ve¬ 
locity  perturbations  in  the  lithosphere  below  the  arrays.  Since  teleseismic  P  is  an  isolated 
(in  time)  arrival  that  travels  near  vertically  through  the  lithosphere  the  effects  of  horizon¬ 
tally  layered  structure  and  scattering  are  effectively  decoupled  and  analytic  infinite  random 
medium  theories  may  be  used.  These  theories  predict  that  forward  scattering  will  predomi¬ 
nate.  that  it  is  controlled  by  the  velocity  fluctuations  (Wu  and  Aki,  1985),  and  that  acoustic- 
theory  may  be  used.  The  measured  quantities  arc’  statistical  variable's  of  the  wavefield  such 
as  variances  of  the  log  amplitude  and  phase,  and  spatial  and  angular  correlations  and  co- 
variances  of  these  quantities.  In  the  most  recent  work  Flatte  and  Wu  (1988)  treat  13  Hz 
waves  at  the  NORSAR  array  (diameter  110  km,  minimum  element  spacing  3  km).  They 


find  the  data  can  be  explained  by  two  overlapping  layers  with  power  law  (i.e.,  scale  free) 
spectra  of  velocity  fluctuations.  The  upper  layer,  from  0--200  km,  has  a  flat  power  spectrum, 
while  the  lower  from  15-250  km  has  a  power  spectrum  that  varies  as  k~‘l  (‘‘self-similar’  i. 
where  k  is  the  wavenumber  of  the  spectrum  of  the  fluctuations.  The  RMS  size  of  the  velocity 
fluctuations  is  a  few  %,  similar  to  other  investigations. 

One  result  of  some  importance  from  this  work  is  the  short  distance  (~  10  km)  over 
which  arrivals  lose  coherence.  For  regional  phases  at  higher  frequencies  even  more  rapid 
loss  of  coherence  might  be  expected.  This  indicates  that  information  from  dense  arrays  of 
small  aperture  is  needed,  and  since  travel  time  and  amplitude  fluctuations  will  be  smaller 
new  analysis  methods  are  desirable.  Accordingly,  to  evaluate  the  effect  of  scattering  on 
direct  phases  the  spatial  coherency  has  been  calculated  using  methods  developed  for  analysis 
of  strong  motion  array  records  (Vanmarcke,  1983).  The  coherency  is  a  frequency  domain 
equivalent  of  the  correlation  function  used  by  Bungum  et  al.  (1985)  and  Ingate  et  al.  (1985). 
Interpretation  of  the  coherency  for  regional  phases  presents  a  challange,  however,  because  of 
the  horizontal  propagation  of  energy  in  a  horizontally  layered  structure  (the  crust):  under 
such  circumstances  the  effects  of  scattering  and  propagation  are  not  decoupled  (e.g.,  Kennett. 
1986).  Two  methods  are  used  to  help  solve  this  problem.  One  is  a  purely  statistical  method 
of  constructing  a  more  or  less  “random”  set  of  seismograms  by  starting  with  an  observed 
trace  and  producing  new  traces  by  randomly  perturbing  the  Fourier  amplitude  and/or  phase. 
While  this  does  not  give  any  information  about  the  earth,  it  allows  us  to  see  what  features 
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of  the  wavefield  are  atfecting  the  coherency.  The  other  method  is  to  compute  theoretical 
seismograms  for  layered  scattering  media  by  finite  differences,  calculate  the  coherency  and 
compare  it  with  observations.  However,  due  to  computational  limitations  it  is  not  possible 
to  make  a  full  comparison. 


METHODS— DATA  ANALYSIS  AND  THEORY 


Computing  Spatial  Coherency 


The  lagged  coherency  is 


C(x,u?) 


Sij{x,ul) 


where  is  the  cross  spectrum  of  seismograms  s(xi),  s(xj)<  lagged  to  remove  time  shifts  due 
to  travelling  wave  propagation;  the  separation  x  —  x,  —  xr  S,,  and  Sjj  are  the  autospectra  of 
the  two  seismograms.  To  compute  the  coherency  we  used  methods  due  to  Jenkins  and  Watts 
(1969)  and  Vanmarcke  (1983).  These  authors  demonstrate  that  gross  time  shifts  due  to  wave 
propagation  across  the  array  must  be  removed  before  accurate  estimates  of  the  coherency 
can  be  made.  Initial  estimates  of  the  lag  times  were  found  by  using  frequency-wavenumber 
analysis  (Capon,  1969)  to  estimate  the  velocity  and  azimuth  of  the  mean  plane  wave  and 
then  calculating  the  lags.  The  cross  correlation  and  autocorrelations  of  a  shifted  pair  of 
seismograms  are  then  windowed  with  a  Hamming  window  of  width  0.5  1  s  and  the  cross 
spectrum,  autospectra  and  coherency  computed.  I  he  purpose  of  the  Hamming  windowing 
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was  to  provide  frequency  smoothing  to  control  bias  at  low  coherencies  (Jenkins  and  Watts, 
1969).  To  see  the  necessity  for  this  note  that  if  no  frequency  smoothing  is  applied,  it  may  be 
shown  that  C  in  (1)  is  identically  unity,  i.e.,  two  perfect  harmonic  components  from  a  finite 
length  signal  are  perfectly  coherent.  This  bias  becomes  more  troublesome  as  the  coherency 
decreases.  We  have  demanded  that  the  frequency-bandwidth  product  be  at  least  10;  this 
allows  successful  calculation  of  coherencies  of  magnitude  0.4-0. 5  and  higher. 

Lengths  of  data  in  the  original  seismograms  are  3-25  s.  This  procedure  is  then  repeated 
for  every  possible  pair  of  seismograms  in  the  array.  In  the  plots  of  coherency  (i.e.,  the 
absolute  value  of  the  coherency)  and  phase  lag  (  he  value  of  the  phase  of  the  calculated 
coherency,  a  complex  quantity),  each  cross  is  such  a  determination.  These  determinations 
are  then  averaged  over  an  appropriate  spatial  interval  and  the  standard  deviations  found, 
using  the  Fisher  Z-transform  (Jenkins  and  Watts,  1969)  for  the  coherency;  these  averages 
are  plotted  as  points  with  error  bars.  Finally,  the  phase  lag  is  checked  to  ensure  that  the 
averaged  values  are  close  to  zero  (within  ±7t/2);  if  they  are  not  a  new  value  of  the  velocity 
and/or  azimuth  is  tried  until  they  are.  We  do  not  adjust  the  lags  other  than  through  the 
velocity  and  azimuth — on  occasion,  there  is  residual  structure  in  the  phase  lags  indicating 
that  the  signal  is  not  a  simple  plane  wave. 
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Finite  Difference  Modelling  in  Highly  Heterogenous  Media 

To  help  analyse  the  results  in  terms  of  scattering  of  seismic  waves  in  media  with  randomly 
varying  velocity  fields,  we  use  the  finite  difference  method  for  producing  synthetic  seismo¬ 
grams.  Coherency  can  then  be  calculated  for  a  set  of  synthetic  seismograms  from  a  model 
with  specified  statistics  and  compared  with  the  observations.  The  technique  is  attractive 
because  it  produces  a  full  solution  to  the  elastic  wave  equation,  thus  all  direct,  converted, 
diffracted  and  guided  waves  are  accurately  modelled.  Also,  unlike  various  high  frequency 
approximations,  there  is  no  limitation  on  the  ratio  of  scatterer  size  to  wavelength.  This 
property  is  essential  when  modelling  waves  in  media  with  self-similar  velocity  fields. 

The  equations  of  motion  for  wave  propagation  in  a  two-dimensional,  elastic,  isotropic 
media  can  be  written  in  terms  of  horizontal  (u)  and  vertical  ( w )  particle  displacements. 


d2u  d 

/,  ^  sdu  ^dw 

d 

du  dw 

P~dW~  di 

{x  +  2ll)—  +  \— 

^dz  +  dx\ 
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d 
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, .  „  ,  dw  .  du 
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In  this  formulation  the  Lame  constants  A  and  p  and  density  p  are  freely  varying  functions 
of  position  ( x,z ),  and  F  is  the  body  force  term.  Unlike  earlier  finite  difference  studies  in 
this  area  (e.g.,  Frankel  and  Clayton,  1986),  we  use  a  staggered  second  order  finite  difference 
scheme  to  solve  Equation  (2).  This  method  is  capable  of  more  accurately  modelling  rapidly 
varying  media,  in  part  because  the  effective  sampling  rate  of  the  velocity  field  is  twice  as 
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great  as  that  for  a  non-staggered  grid  (Virieux,  1984,  1986).  At  present  we  have  limited  our 
computations  to  two-dimensional  models. 


Constructing  Random  Velocity  Perturbations 


Random  media  can  be  characterized  by  three  basic  quantities,  their  correlation  function 
(or  fluctuation  spectrum),  their  probability  distribution,  and  their  standard  deviation  (or 
spectra]  amplitude).  There  are  many  different  correlation  functions  which  can  be  used  to 
describe  a  random  medium,  but  in  this  presentation  we  will  only  be  concerned  with  zeroth 
order  Von  Karman  functions  (Tatarskii,  1961).  In  two  dimensions,  the  wavenumber  domain 
representation  of  this  class  of  functions  can  be  written  in  terms  of  the  correlation  length  a 
and  the  horizontal  and  vertical  wavenumbers  kx  and  kz  as, 


P(kM 


(3) 


The  outstanding  feature  of  media  with  this  correlation  function  is  the  presence  of  roughness 
at  all  length  scales.  This  property  is  often  seen  in  nature  (Mandelbrot,  1977),  and  is  thought 
to  best  reflect  the  random  portion  of  the  earth’s  velocity  field  (Wu,  1982;  Dainty,  1984; 
Frankel  and  Clayton,  1986). 

Two-dimensional  realizations  for  each  of  the  three  material  properties,  Lame’s  parameter, 
shear  modulus,  and  density,  were  constructed  by  inverse  Fourier  transforming  the  product 
of  the  wavenumber  representation  of  the  zeroth  order  Von  Karman  function  and  a  random 
phase  term.  When  transformed  back  to  the  space  domain,  each  of  the  three  functions  will 
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have  the  desired  autocorrelation  function,  zero  mean  and  a  Gaussian  probability  distribution. 


OBSERVATIONS  AT  FENNOSCANDIAN  ARRAYS  FROM 

QUARRY  BLASTS 

To  examine  the  coherency  of  Lg  from  regional  events  data  from  the  three  Fennoscandian 
arrays  mentioned  previously,  ARCESS,  FINESA  and  NORESS,  was  used.  The  locations  of 
the  arrays  are  given  in  Table  1  and  Figure  1  shows  the  configurations.  While  all  three  arrays 
contain  three  component  stations,  only  vertical  components  have  been  used  in  this  study. 
Further  details  of  the  arrays  are  in  the  cited  references.  Results  are  presented  for  four  events 
in  the  distance  range  200-400  km;  information  about  these  events  is  summarized  in  Table  2. 
All  of  the  events  are  believed  to  be  quarry  blasts  either  on  the  basis  of  location  at  a  known 
quarry  and/or  the  presence  of  a  strong  Rg  phase.  Time  windows  for  analysis  are  selected  to 
include  not  only  Lg  but  also  Sn  and  Sg,  if  they  are  present.  Any  Rg  phase,  however,  was 
excluded.  It  was  found  that  the  coherency  analysis  was  not  sensitive  to  small  shifts  of  the 
window. 

To  illustrate  the  analysis,  we  choose  the  15  November  1985  event  recorded  at  the  FINESA 
array.  Figure  2  shows  the  array  seismograms  with  the  analysis  window  indicated.  In  Figure 
3  the  power  spectrum  and  the  frequency-wavenumber  diagram  for  a  frequency  of  3.6  Hz  are 
presented.  Figures  4  and  5  show  the  absolute  value  of  the  coherency  and  the  phase  of  the 
coherency  of  the  lagged  traces.  Figure  5  demonstrates  that  the  averaged  phase  lags  are  within 
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±tt/2  of  zero,  as  required.  Figure  4  illustrates  that  coherency  decays  with  increasing  spatial 
separation  more  rapidly  at  higher  frequencies.  Figure  6  shows  the  seismograms  at  the  array 
center  seismometer  for  the  other  events,  and  Figure  7  presents  the  absolute  value  of  coherency 
for  the  event  recorded  at  NORESS,  illustrating  more  clearly  the  increasingly  rapid  decline 
of  coherency  w'ith  separation  as  frequency  increases.  In  Figures  4,  5  and  7,  and  throughout 
this  paper,  spatial  separation  is  taken  as  the  absolute  distance  between  seismometer  pairs 
without  regard  to  whether  the  separation  vector  is  oriented  along  or  perpendicular  to  the 
azimuth  of  the  wave.  We  find  that  there  is  no  apparent  difference  between  these  cases,  in 
agreement  with  w'ork  on  teleseismic  arrivals  (Aki,  1973;  Capon,  1974)  and  strong  motion 
(Harichandran  and  Vanmarcke,  1986).  This  is  contrary  to  the  findings  of  Der  et  al.  (1984). 
but  most  of  these  authors’  results  are  at  lower  frequencies. 

As  stated  above,  coherency  declines  with  increasing  separation  and  declines  faster  for 
higher  frequencies.  We  attribute  this  decline  in  coherency  to  scattering.  By  looking  at 
different  frequencies  we  have  examined  different  wavelengths.  In  an  attempt  to  gain  an 
integrated  look  at  the  data,  Figures  8  and  9  plot  the  spatially  averaged  coherencies  for  all 
events  and  frequencies  analysed  against  distance  in  wavelengths.  Intriguingly,  within  error 
limits  it  appears  that  the  curves  of  decline  of  coherency  with  separation  for  the  different 
frequencies  are  all  very  similar  if  the  separation  is  scaled  to  the  wavelength.  The  effect  is 
most  dramatic  for  the  event  recorded  at  NORESS.  The  curves  for  different  events  appear 
to  be  similar;  this  may  reflect  the  similarity  of  site/path  heterogeneities  for  these  crystalline 
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rock  sites  within  the  Fennoscandian  shield.  Note  that  the  coherency  declines  to  ~  0.5  within 
about  one  wavelength. 

For  a  preliminary  interpretation  of  these  results,  we  note  that  in  spite  of  travelling  several 
hundred  kilometers  through  the  lithosphere  regional  seismograms  still  show  distinct  phases 
such  as  P  and  Lg,  i.e.,  the  effect  of  scattering  is  not  so  strong  as  to  obliterate  all  features 
of  wave  propagation.  This  suggests  the  scaling  illustrated  in  Figures  8  and  9  must  be  a 
propagation  effect  and  indicates  that  the  heterogenities  in  the  lithosphere  must  also  scale 
in  some  manner.  This  has  led  us  to  investigate  self-similar  models  of  velocity  and  density 
perturbations  in  the  crust  as  discussed  in  the  next  section.  The  term  “self-similar”  is  here 
taken  to  mean  having  a  power  law  spectrum,  at  least  within  a  certain  range  of  length  scales. 
For  strict  self-similarity  this  range  is  0  to  oo,  but  it  is  not  possible  to  construct  a  finite 
realization  for  such  a  case  because  a  power  law  spectrum  will  tend  to  oo  at  one  of  the  limits. 
Thus  the  Von  Karman  function  defined  in  Equation  (3)  is  asymptotically  self-similar  for 
short  length  scales  (high  wavenumbers)  but  is  not  self-similar  for  length  scales  longer  than 
the  parameter  a,  which  is  a  length.  The  appropriate  length  scales  are  set  by  the  wavelengths 
of  the  waves  analysed,  which  are  in  the  range  0.4-2  km. 
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THEORETICAL  INVESTIGATIONS 


Coherency  of  Randomized  Time  Series 

In  interpreting  coherency  it  is  of  interest  to  know  whether  fluctuations  in  amplitude  or 
phase  (travel  time)  of  the  seismograms  are  important.  To  examine  this  question  randomized 
time  series  were  constructed,  starting  with  real  data,  and  then  analysed  as  if  they  were 
spatially  separated  seismograms.  An  example  is  shown  in  Figure  10.  The  top  trace  is  the 
top  trace  of  Figure  2  within  the  window  shown  on  Figure  2.  To  produce  the  second  trace  the 
Fourier  transform  of  the  top  trace  is  calculated,  randomly  perturbed  while  maintaining  its 
symmetry  properties  about  zero  frequency,  and  transformed  back  into  time.  In  this  case  only 
the  phase  of  the  Fourier  transform  was  perturbed;  Figure  11  shows  the  perturbations,  which 
are  scaled  with  frequency  so  as  to  model  travel  time  shifts  of  0.04  s  maximum.  The  third 
trace  is  constructed  from  the  second  trace  in  a  similar  manner  (with  a  different  realization 
of  perturbations)  and  so  on,  i.e.,  a  set  of  “seismograms”  that  are  succesively  more  unlike 
the  initial  seed  are  produced,  just  as  might  be  expected  for  a  suite  of  records  at  increasing 
spatial  separation.  Figure  12  presents  coherency  at  4  Hz  for  the  time  series  shown  in  Figure 
10  computed  as  if  successive  traces  were  200  m  apart.  The  coherency  decays  with  increasing 
“separation”  and  is  of  the  same  order  of  magnitude  as  the  observations  presented  earlier. 
Numerical  experiments  such  as  this  have  demonstrated  that  phase  perturbations  are  much 
more  important  than  amplitude  perturbations  in  controlling  coherency. 


Finite  Difference  Modelling 


While  the  effort  using  randomized  fields  is  instructive,  it  is  also  clearly  not  reproducing 
some  features  of  the  data.  To  give  one  example,  the  seismograms  in  Figure  2  all  exhibit 
the  characteristic  of  a  definite  “S”  arrival  at  about  29  s.  In  Figure  10,  this  corresponds  to 
the  peak  at  1.8  s  on  the  first  trace,  which  is  taken  from  the  top  trace  of  Figure  2,  but  the 
other  traces  do  not  in  general  maintain  this  characteristic.  The  fundamental  difficulty  is 
that  the  “model”  is  not  physical  and  thus  does  not  contain  the  concept  of  seismic  phases, 
for  example.  To  remedy  this  we  have  begun  to  use  finite  difference  calculations  to  compute 
array  seismograms  which  can  then  be  analysed  in  the  same  manner  as  observed  seismograms. 
However,  while  the  finite  difference  method  allows  the  calculation  of  complete  seismograms 
in  a  physical  model,  computational  limitations  cause  the  models  that  can  be  realistically 
considered  to  be  smaller  than  desirable  and  the  signals  to  have  a  narrow  frequency  range. 

The  work  in  finite  difference  lias  just  begun  and  so  only  an  example  will  be  presented 
here.  In  Figure  13  a  two-dimensional  mwdel  of  a  layer  over  a  half  space  with  a  rough  interface 
and  random  velocity  and  density  perturbations  is  shown.  For  these  trial  calculations  we  have 
limited  the  size  of  the  model  (6.5  x25  km)  and  the  frequency  content  of  the  seismic  signal  (1- 
5  Hz).  The  average  model  has  a  6  km  thick  layer  of  P  velocity  6  km/s  overlying  a  substrate 
with  P  wave  velocity  7  km/s.  Poisson’s  ratio  is  0.25.  A  power  spectrum  for  the  fluctuations 
of  the  form  (3)  was  used  with  the  parameter  a  =  1.5  km,  equal  to  the  wavelength  of  P 
waves  at  the  mean  frequency  of  3  Hz.  The  RMS  velocity  fluctuations  are  0.15.  Synthetic 
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seismograms  for  the  vertical  component  arc  shown  in  Figure  13  with  an  analysis  window  for 
P  indicated.  P  was  chosen  for  analysis  because  the  source  is  an  explosion  and  the  distance 
travelled  is  relatively  short,  about  10  wavelengths.  These  seismograms  were  taken  for  an 
array  between  22  and  24  km  distance  on  the  surface  of  the  model  and  were  spaced  125  m 
apart.  The  appropriate  lagging  velocity  for  these  seismograms  in  this  window  was  found  to 
be  4.5  km/s,  slower  than  the  P  wave  velocity  of  6  km/s.  This  indicates  scattering  to  S. 

Figure  14  shows  the  coherency  at  the  center  frequency  of  3  Hz  for  these  synthetic  seismo¬ 
grams.  Noting  that  the  wavelength  at  4.5  km/s  and  3  Hz  is  1.5  km,  we  see  that  the  results 
in  Figure  14  agree  roughly  with  Figures  4  and  7.  A  further  comparison  emphasising  the 
agreement  is  shown  in  Figure  15  which  presents  the  synthetic  results  on  a  wavelength  scaled 
plot  with  the  observations  from  the  NORESS  event.  Work  with  layered  media  for  strong 
motion  modelling  (Toksoz  et  al.,  1989)  has  demonstrated  that  the  velocity  perturbations  are 
more  important  than  the  rough  interface  in  producing  a  loss  of  coherency. 

DISCUSSION  AND  CONCLUSIONS 

A  crucial  issue  in  research  into  seismic  scattering  is  the  nature  of  the  earth  structure  causing 
it.  In  this  investigation  we  have  concentrated  on  the  seismic  S  phases  at  regional  distances; 
this  limits  the  investigation  to  the  crust  (Lg)  or  possibly  the  uppermost  mantle  (Sn).  By 
choosing  locations  in  the  crystalline  province  of  Fennoscandia  any  effects  due  to  near-surface 
sedimentary  layers  has  been  eliminated.  Also,  by  looking  at  direct  phases  we  have  empha- 
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sized  forward  scattering  as  opposed  to  backscattering.  Forward  scattering  is  sensitive  to 
large  scale  fluctuations  and  to  velocity  perturbations.  Backscattering  is  sensitive  to  smaller 
scale  fluctuations  than  forward  scattering  for  the  same  wavelength  of  the  incident  wave  and 
to  impedance  contrasts  rather  than  velocities  (Wu  and  Aki,  1985).  Thus  the  scatterers  re¬ 
sponsible  for  the  effects  seen  in  this  study  are  not  necessarily  the  same  as  those  responsible 
for  such  phenomena  as  coda  and  attenuation  which  are  controlled  by  backscattering. 

Our  array  measurements  of  coherency  in  Fennoscandia  demonstrate  that  it  declines  to 
~  0.5  within  about  a  wavelength.  This  wavelength  scaling  is  scaling  according  to  the  length 
scale  of  the  wave;  there  is  no  sign  in  the  wavelength  range  examined  (0.4-2  km)  of  a  length 
scale  for  the  medium  afFecting  the  results.  This  is  explainable  by  propagation  in  a  medium 
with  a  power  law  spectrum  of  medium  heterogeneity  such  as  those  used  by  Flatte  and  Wu 
(1988)  for  teleseismic  observations.  The  coherency  statistics  appear  to  be  isotropic  (i.e.,  they 
do  not  depend  on  direction  with  respect  to  the  azimuth  of  the  incoming  wave),  in  agreement 
with  observations  of  teleseismic  phases  in  the  same  region  (Flatte  and  Wu,  1988)  and  strong 
motion  observations  in  Taiwan  (Harichandran  and  Vanmarcke,  1986).  We  do  not  observe 
the  distinction  between  on  azimuth  and  perpendicular  to  azimuth  separations  seen  at  lower 
frequencies  by  Der  et  al.  (1984). 

To  assess  the  causes  of  the  loss  of  coherency,  randomized  time  series  were  examined  to 
see  what  features  of  the  observed  wavefield  were  important.  Phase  (or  travel  time)  variation 
was  found  to  be  dominant  over  amplitude  variation,  reinforcing  the  importance  of  velocity 
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perturbations,  since  they  control  travel  times.  Travel  time  variations  of  the  order  of  0.01  Rec¬ 
over  100  m  at  4  Hz  were  found  to  be  sufficient  to  explain  the  data. 

To  provide  a  satisfactory  explanation  of  the  observations,  a  physical  model  of  the  prop¬ 
agation  of  seismic  waves  through  appropriate  earth  models  which  include  scatterers  is  re¬ 
quired.  In  preliminary  work,  similar  efFects  to  those  observed  have  been  found  in  synthetic 
seismograms  calculated  for  models  with  rough  boundaries  and  random  variations  of  veloc¬ 
ity  and  density.  At  present,  however,  computational  limitations  have  prevented  a  thorough 
test.  The  linear  dimensions  of  the  models  are  too  small  for  direct  comparison  and  they  are 
two-dimensional  rather  than  three-dimensional.  Related  to  the  problem  of  model  size  is  the 
high  value  of  RMS  perturbation  that  had  to  be  used,  15%,  to  produce  the  desired  loss  of 
coherency.  This  is  considerably  higher  than  the  1-5%  range  found  in  studies  of  teleseismic 
data  in  the  same  region  (Capon  and  Berteussen,  1974;  Flatte  and  Wu,  1988),  although  these 
studies  are  at  somewhat  lower  frequencies.  In  addition,  we  have  not  yet  tested  the  wave¬ 
length  dependence  of  the  loss  of  coherency.  Future  work  will  concentrate  on  this  question 
both  by  using  more  powerful  computers  for  the  finite  difference  calculations  and  by  exploring 
analytical  methods  such  as  those  of  Flatte  et  al.  (1979). 
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Table  1:  Array  locations. 


Array 

Latitude 

Longitude 

No.  of  Elements 

ARC  ESS 

69.535N 

25.506E 

25 

FINESA 

61.444N 

26.079E 

10 

NORESS 

60.735N 

11.541E 

25 

Table  2:  Events  Analysed.  Azimuth  (from  N)  and  Velocity  are  estimated  from 
wavenumber  analysis,  not  locations. 

frequency- 

Date 

OT 

(UTC) 

L&t. 

Long. 

Location 

Recorded 

Dist. 

(km) 

Az. 

Velocity 

(km/s) 

27  Feb.  1985 

12:58:31 

59.3N 

6.5E 

Blasj0“ 

NORESS 

324 

243° 

3.8 

15  Nov.  1985 

13:53:37 

61. IN 

29.9E 

Leningrad6 

FINESA 

209 

99° 

3.6 

6  Feb.  1986 

12:22:04 

59.3N 

28. IE 

Leningradc 

FINESA 

264 

151° 

4.5 

16  Dec.  1987 

11:48:54 

68. IN 

33.2E 

Kolad 

ARCESS 

350 

116° 

4.1 

flBla6j0  Quarry  location  from  Dysart  and  Pulli  (1987).  OT  from  S — P  time. 
JOT  and  location  from  Ahjos  ei  at.  (1986a) 
cOT  and  location  from  Ahjos  ti  at.  (1986b) 

dQuarry  location  from  Mykkeltveit  et  a!.  (1987).  OT  from  S- — P  time. 
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FIGURE  CAPTIONS 


Figure  1:  Geophone  array  configurations  at  FINESA  (top)  and  NORESS  and  ARCESS 
(bottom) 

Figure  2:  Vertical  component  seismograms  recorded  at  FINESA  15  Nov.  1985  from  a  quarry 
blast  at  209  km  range. 

Figure  3:  The  squared  Fourier  spectrum  and  the  frequency-wavenumber  analysis  for  the 
windowed  portion  of  the  data  in  Figure  2. 

Figure  4:  Coherency  vs  distance  at  2  Hz  (top),  4  Hz  (middle)  and  6  Hz  (bottom)  for  the 
windowed  data  in  Figure  2. 

Figure  5:  The  phase  lag  vs  distance  at  2  Hz  (top),  4  Hz  (middle)  and  6  Hz  (bottom)  for 
the  windowed  data  in  Figure  2. 

Figure  6:  Vertical  component  seismograms  recorded  at  FINESA  6  Feb.  1986  (top),  ARCESS 
16  Dec.  1987  (middle)  and  NORESS  27  Feb.  1985  (bottom)  center  seismometers  with  analysis 
windows  indicated. 

Figure  7:  Coherency  vs  distance  at  2  Hz  (top),  4  Hz  (middle)  and  6  Hz  (bottom)  for  the 
NORESS  windowed  data  in  Figure  6. 

Figure  8:  Coherency  vs  distance  normalized  to  the  wavelength  for  events  recorded  at 
FINESA. 

Figure  9:  Coherency  vs  normalized  distance  for  events  recorded  at  ARCESS  (top)  and 
NORESS  (bottom). 
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Figure  10:  Randomized  time  series.  Top  trace  is  taken  from  analysis  window  of  top  trac^ 
of  Figure  2.  Succeeding  traces  are  produce  by  iterative,  random  perturbations  of  Fourier 
phase  of  top  trace.  See  text  for  discussion. 

Figure  11:  Fourier  phase  perturbations  (radians)  used  to  produce  second  trace  from  the 
top  in  Figure  10  from  the  top  trace.  Perturbations  are  added  to  the  phase.  See  text  for 
discussion. 

Figure  12:  Coherency  vs  distance  for  the  randomized  traces  in  Figure  10  at  4  Hz.  The 
distance  scale  was  constructed  by  considering  each  trace  in  Figure  10  to  be  part  of  a  linear 
array  of  spacing  200  m. 

Figure  13:  (Top)  Model  used  for  synthetic  computations.  The  source  is  explosive  and 
seismograms  are  computed  for  an  array  2  km  long  with  a  spacing  of  125  m.  Dimensions  of 
the  model  are  6.5  x  25  km;  the  upper  layer  has  a  P  velocity  of  6  km/s,  the  lower  7  km/s  and 
Poisson’s  ratio  is  0.25.  (Bottom)  Vertical  seismograms  for  a  wavelet  of  dominant  frequency 
3  Hz,  lagged  at  4.5  km/s  and  with  an  analysis  window  indicated.  Top  number  to  right  of 
each  trace  is  the  element  number  (01  is  at  22  km,  17  at  24  km)  and  the  bottom  number  is 
the  maximum  ampllitude. 

Figure  14:  Coherency  vs  distanceat  3  Hz  for  the  vertical  component  synthetic  seismograms 
of  Figure  13. 

Figure  15:  Coherency  vs  normalized  distance  for  the  27  Feb.  1985  event  recorded  at 
NORESS  and  the  synthetic  seismograms  of  Figure  13.  See  text  for  discussion. 
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Abstract 


The  existence  of  Rayleigh-scattering  of  elastic  plane- waves  by  anisotropic  homoge¬ 
neous  inclusions  is  theoretically  demonstrated.  The  case  of  transverse  isotropy  is  stud¬ 
ied  in  detail.  It  is  shown  that,  a  scattered  longitudinal  wave  creates  radial-longitudinal 
(P),  collatitudinal-shear  (SV)  and  azimuthal  shear  waves  (SH).  Likewise,  scattered  SV 
and  SH  waves  generate  each,  radial  P  waves,  SV  and  SH  shear  waves.  All  scattered 
amplitudes  are  proportional  to  the  square  of  the  frequency  and  have  radiation  pattern 
signatures  as  those  of  equivalent  dipoles,  center  of  compression  and  double  couples.  It 
is  shown  that  observations  of  spatial  patterns  of  scattered  amplitudes  can  yield  through 
inversion,  the  elastic  constants  of  the  anisotropic  inclusion.  Our  results  can  serve  as 
a  theoretical  basis  for  the  observed  short-period  SH  and  SV  waves  from  underground 
explosions  at  teleseismic  distances. 
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INTRODUCTION 


Elastic  wave  scattering  by  elastic  heterogeneities,  using  the  equivalent  source  method  and 
Born  approximation  has  been  successfully  formulated  and  applied  to  acoustics  and  seismology 
(Gubernatis  et  al.,  1977;  Haddon  and  Cleary,  1974;  Wu  and  Aki,  1985).  Recently,  Ben- 
Menahem  (1989a,  1989b)  obtained  the  explicit  analytic  Green’s  tensor  with  the  ensuing 
displacements  fields  due  to  unipolar  and  dipolar  point  sources  in  azimuthally-isotropic  media. 
In  the  present  paper,  we  combine  the  results  of  the  above  two  disciplines  to  show  that  the 
method  can  be  generalized  to  include  scattering  due  to  a  combined  elastic  heterogeneity 
and  anisotropy  having  slightly  different  elastic  constants  from  the  surrounding  medium.  In 
particular,  we  work  out  in  detail  the  scattering  due  to  a  small  azimuthally-isotropic  inclusion. 
We  show  that  scattered  field  is  equivalent  to  those  generated  by  point  dipoles  and  double 
couples.  The  most  important  results  of  our  study  is  that  a  few  scattering  experiments, 
through  which  plane  P,  SH  and  SV  waves  are  sent  to  interact  with  the  inclusion,  lead  to 
the  determination  of  the  elastic  constants  of  the  inclusion.  We  hope  that  this  will  become 
useful  in  exploration  geophysics  where  azimuthal-isotropy  is  linked  to  the  presence  of  shale 
(Winterstein,  1986).  Also,  it  may  explain  the  presence  of  SH  and  SV  waves  in  the  near  and 
far  fields  of  underground  chemical  and  nuclear  explosions  (Kisslinger  et  al.,  1961;  Gupta  and 
Blandford,  1983). 


1.  EQUIVALENT  SOURCES  FOR  FIRST-ORDER 


SCATTERERS 


Consider  an  homogeneous  anisotropic  elastic  solid  of  density  p.  The  force- free  Cauchy  equa¬ 
tion  of  motion  for  this  solid  is  (e.g.,  Ben-Menahem  and  Singh,  1981), 

-»  d2u 

dlvT~pdF  =  0  (L1) 

where  T  ( u )  is  the  stress  tensor,  and  u(r,  t)  are  the  spectral  displacements. 


We  choose  a  Cartesian  coordinate  system  0(x,  y,  z)  in  which  xy  is  the  plane  of  isotropy. 
In  matrix  notation,  the  stress  tensor  for  azimuthally  isotropic  media  bears  the  explicit  matrix 
form  (Love,  1944) 


(1.2) 


Here,  (ux,  uy,  uz)  are  the  Cartesian  displacement  components  and  { p ,  C u,  CJ3,  C33.  C44, 
C66j  C\2  =  Ci\  —  2Cee}  are  six  independent  structural  parameters  of  the  medium.  The 
matrix  relation  in  Eq.  (1.2)  can  be  put  in  the  convenient  tensor  form 


T=  [C, 3  7  div  u  +  C,4(Vu  +  uV)]+  T„ni«,tT0pu- 


(1.3) 
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(dxi  d a  \ 

~dx  ~dy  J 

du 

+  (C33  —  C13  —  2C44)  ezez 


(1.4) 


^dux  _  _  nd-uv  _  _  ( dux  duy\  _  _  _  . 

2  +  ~d^exCx  -  {'W  +  ~fa)  {ex6y  +  eyei) 


+  (C44  —  Ce,&) 

Consider  an  anisotropic  obstacle  of  structural  parameters  {/);  Cn]  Ci3;  C33;  C44;  C66} 
embedded  in  an  infinite  isotropic  elastic  space  with  parameters 


C[?  =  cg>  =  A  +  2/i  ;  Cff  =  eg  =  ft  ;  (C13)°  =  A  . 


Let 


C,7  =  40)  + 


(1.5) 


where  6C7,  are  small  perturbations  such  that  C{j  =  C,-°'  in  the  isotropic  limit.  We  assume 
that  the  obstacle  is  small  compared  with  the  wavelength  of  the  impinging  radiation.  We  put 
Eqs.  (1.3)  and  (1.4)  into  Eq.  (1.1)  with  CtJ  given  by  Eq.  (1.5),  and  assume  a  perturbation 
solution  of  Eq.  (1.1)  in  the  form  u  =  u ^  4-  u1,  where  is  a  primary  field  and  is  a 
scattered  field.  The  equation  of  motion  then  splits  into  two  simultaneous  wave  equations, 

div  {C\l]I  div  u(0)  +  C[?{Vd0)  +  u<0)V)}  -  =  0  0-6) 

div  {Cl?I  div  u<‘>  +  CfffVtZW  +  tZ^V)}  -  P =  ~div  M  (1.7) 

oz 


where 


M=  [bCu  —  SC\3  —  26C44] 


45 


gu(o)  ( <9u(0>  du(0)\ 

+  [6C33  —  6C\z  —  2<5C44]  —q^ - [«Cn  —  SC\3  —  6C44]  (  —— — I — -- J  e2e2 

du<°)  /5u(°)  du<0>\  ^ 

+  [6C44  —  <5C66]  2 eyey  +  2—^—exex  —  I  — ■— — I — ~faTj  "h  ^ve*  + 

+  <5C44[Vu<0)  +  u(0)V]  .  (1.8) 

The  right-hand  side  of  Eq.  (1.7)  indicates  the  presence  of  a  dipolar  point  source  with  moment 
density  M  (Ben-Menahem  and  Singh,  1981,  pp.  167-169;  Tarantola  et  ah,  1988).  Through 
this  representation,  M  is  identified  as  the  equivalent  source  of  the  scattered  field.  Indeed, 
introducing  the  explicit  expression  for  [Vu*0'  +  u^V],  we  may  decompose  M  in  Eq.  (1.8) 
into  five  partial  tensors,  each  of  which  represents  a  different  force  system: 

I.  Explosion 

div  u W  0  0 

M,=  (SCI3)  q  dw  ~(0)  Q  (1.9) 

0  0  div  u(0) 


46 


II.  Dipole  along  the  symmetry  axis  (2) 


0  0  0 


A/2—  (^(^33  ~  <^Cl3) 


0  0 


0 


0  0  dzu (°) 


III.  Dipoles  in  the  isotropy  plane  ( xy ) 


M  3—  (6Cn  —  8C\  3) 


dxug>  +  dvuW  0  0 


dxu<?)  +  dyUw  0 


dyvtp  0  0 


-  2 (tfCec) 


0  dxuW  0 


0  0  0 


(1.10) 


(1.11) 
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IV.  Two  double-couples  on  perpendicular  saggital  planes  (xz,  yz) 


dxu[0)  +  8zu{0) 


M  4—  (^^44) 


a,u[° >  +  d,uW  (U2) 


dxu(^  +  dzu(°]  dyu[°}  +  dzu^)]  0 


V.  Double-couple  in  the  isotropy  plane  ( xy ) 


dtuW+dyUM  0 


a/5=(^'6C)  a,«s?)  +  dA* 


(M3) 


In  the  isotropic  limit 


<5C-i3  —  6\  ;  8 C 33  —  8C 13  —  8Cu  —  8C 13  —  28  ;  8C  4 4  —  8C gg  —  8y.  ,  (1.14) 


and  we  fall  back  on  the  known  case  of  scattering  due  to  a  localized  inhomogeneity  (\Vu  and 
Aki.  1985).  To  this  end,  we  can  subdivide  M  in  an  alternative  form,  through  which  this 


48 


limit  is  immediately  apparent.  Using  superscripts,  we  have: 


M 


=  8C 


13 


drUi0)  +  dyltf  +  d,U{0) 


dxUw  +  dyuf  +  dzu{0) 


dtuW  +  dyu™  +  0,u(°> 


(1.15) 


A/'-  6CU 


aru<°)  ^4°>  +  5y4°)  5rui0)  +  ^40) 


^W*01  +  5ywi0)  2(9vuj/0)  <9yn(0)  +  d2?d0) 


dxdu^  +  dzu[0'1  dyuf1  +  duzu[0)  2dzu[0) 


(1.16) 


M  —  (^Cn  —  £C\3  —  2^  C44) 


d'UW+dyUM 


0  0 


$A0)  +  ^  0 


0  0 


(1.17) 
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0  0  0 


M  —  {SC33  —  SC  13  —  2SC44)  Q  q  q 


(1.18) 


0  0  dzuW 


2  dyu[°) 


~(dxu<?)  +  dyUW)  0 


M  =  {SC44  -  SC6e)  _(dxti(o)  +  dyu^)  2dxu<?> 


(1.19) 


~.3  4  ^,5  _  ,-.1  ^2 

In  the  isotropic  limit  A1  ,  M  ,  M  vanish  and  M  >  Af  survive  to  yield  the  scattering  field 
due  to  a  localized  inhomogeneity. 
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2.  THE  SCATTERED  FIELDS 


For  any  given  the  scattered  field  u l1)  solves  Eq.  (1.7)  and  is  given  approximately  by 
Ben-Menahem  and  Singh  (1981) 


(iso) 

U(l)  =  -[M  :  V  G  KM 


(2.1) 


*~(iso) 

where  G  is  the  Green’s  tensor  of  the  isotropic-homogeneous  background  matrix  of  the 
inclusion 


-.(iso) 

G  = 


1 

4  n/j. 


I 


e~ikgR 

R 


+  VV 


e~ikgR  _  e~ikaR' 

M  . 


with 


(2.2) 


,  =  CS>;  A  =  C<?; 


and  ( Sv )  is  the  ‘volume’  of  the  obstacle. 


Since  VV/(i?)  _  (escs  +  e^)^  —  +  eReR^,  in  spherical  coordinates  ( R ,  9,  <f>),  we 
shall  have  in  the  ‘far- field’ 


♦-►(iso) 

47T/J  G  = 


/?2  c~'kaR_  ^  e-ik»R 

^2  '  ^  ^  '(egeg  +  e^e*)  +  0 


1 

k20R 2 


(2.3) 


Substituting  from  Eq.  (1 .9)— ( 1. 13)  into  Eq.  (2.1),  we  obtain  the  five  partial  scattered  fields, 


=  ( SCi3)(div  u(G))div  G  (2.4) 

-<i)  _  ,,n  /<9i40)\  d  ~ 

u2  -  (6C33  ~  SC13)  y~^~j  fc(G  -e2)  (2.5) 
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^1) 

u3 


/duW 

duf\ 

f  9 

-x  d  «  ..1 

1  X 

V  dx 

+  dy  ) 

|aI(G 

6x  +  dy  G  6  J 

= 


41} 


G=G 


♦  (iso) 


(2.6) 


(2.7) 

(2.8) 


(a)  Incident  P-wave  (a  =  yj{\  +  2n)/p) 


Let 


=  Ukp)pe‘“[t~°{ 5'*in<’«+ic08<9«>]  =  ulp  , 


(2.9) 


where  p  is  the  wave’s  normal  and  do  is  the  angle  that  p  makes  with  the  symmetry  axis 
a  (Figure  1).  Due  to  the  azimuthal  symmetry  of  the  scattering  obstacle,  we  can  choose, 
without  loss  of  generality 


p  =  ey  sin  Oo  +  ez  cos  90 


(2.10) 
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Suppressing  the  factor  e,w(,  the  partial  wave  derivatives,  needed  for  Eqs.  (2.4)-(2.8),  are: 


u(_0)  = 


H0) 

dy 

dui0) 

dx 


du ^  du (?)  du^  ~(<*°) 

0  ;  ~~z~~  =  =  0  ;  div  G  = 


dx 

—  iup0kQ  sin2  0o  ; 


3?/  dz 

duM 


-VI- 


— «fca  R ' 


4t(A  +  2/i)  \R 
= —iu^ka  cos2  9q  ;  div  = —iuvQka  (2.11) 


duj0> 

dx 


=  0  ; 


02 

5U(°)  av(°) 

=  -7p-  =  —iuQka  sin  0O  cos  0o  • 
dy  dz 


In  the  far-field ,  we  obtain  from  Eq.  (2.3),  to  terms  of  order  ( kaR ) 


-2 


div  G=  — 


ikn 


.-tk„R 


47r(A  -f  2y)  R 


-cr  . 


(2.12) 


Also, 


_  d 

Cldx 
„  d 
tydy 


[eR  sin2  0  cos2  <f>  +  eg  sin  9  cos  6  cos2  (j>  —  e#  sin  9  cos  <)>  sin  <^>] 


—  =  [e/?sin20sin2<?!>  +  eg  sin  0  cos  0  sin2  0  +  e,*,  sin0cos  <£sin  0} 


d_ 
dR  ' 
_d_ 
dR  ’ 


dz 

-  d  _  d 

er%  +  Cyax 

-  9  _  a 

Cyc>2  +e^t/ 

_  d  _  d 

Cz  r,  ^ X  r\ 

dx  dz 


d  d 

ez  —  =  [eH  cos2  9  —  e#  sin  9  cos  0]  — —  , 
A-  dR 


d 


(e/?  sin2<?isin0  +  eg  sin  2<^cos  0  +  e#  cos2</>]  sin  0  — —  , 

dR 

(eft  sin  29  sin  4>  +  eg  cos  29  sin  <f>  4-  e#  cos  0  cos  <?!>]  r—  , 

C'./t 

[e^sin  20  cos  4>  +  e#  cos  20  cos  4>  —  cos  0  sin  <^7775-  , 

OR 


(2.13) 


where  (R,  0,  (^>)  are  spherical  coordinates  of  the  observation  point  and  (cr,  eg,  e$)  are  the 
associated  unit  vectors  (Fig.  1). 


Using  Eqs.  (2.11)— (2.13),  the  final  expressions  for  the  spherical  displacements  compo- 
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nents,  yielding  the  respective  P,  SV  and  SH  wave  motions,  are 


=  — 06  sin2  sin  0  sin  2<f>  -f  - — —  sin2  0o  cos  0  cos  <p  ,  (2.19) 

fi  // 

are  the  ensuing  radiation-patterns.  Figures  2-6  display  three-dimensional  amplitude  pat¬ 
terns  resulting  from  the  interaction  of  the  obstacle  with  an  incident  P-wave:  FRP(0o  =  0)  in 
Figure  2;  F6p(90  =  0)  in  Figure  3;  FRP(0o  =  90°)  in  Figure  4;  FgP(9Q  =  90°)  in  Figure  5; 
FMOo  =  90°)  in  Figure  6.  Numerical  values  of  the  paramaters  involved  in  the  calculation 


are  given  in  Table  1. 


(b)  Incident  SV-Wave  (d  =  \J p/ p} 


Let 


u<0)  =  ifV){eT  x  p)cM*~ ^(*»in «o+* c°. »o)]  =  u[)sv)(e-x  x  p)  t 


with  the  partial  derivatives 


i 


u'j'  =  0  ; 


du<°)  0U<°>  du<°>  du<°>  dul°> 

dx  ~  dy  ~  dz  ~  dx  ~  dx  ~0, 


•  (sv)i  •  n  ft  ^uy  ■  (SV).  2n 

-r —  =  iu0  ’  kpsmVocosOo  ;  ■—  =  -iKn  ’kn  cos  O0  =  — — 

dy  dy 

du<0) 


=  —  h*oSl  sin 0O  cos  0O  ;  div  u(0'  =  0  . 


dz 


The  corresponding  displacement  components  are: 


rr(SV)kl(t>v)  ft  etu(l 

UR  =  -£Ao  — ^ —  —  ^ft(sv)(0o;  0\  F,  Cn) - — - 

Tr(SV)kl(f>v)  e,u,('  ^ 

ue  =  Uq  — — — Fo(sv){00;  0\  <p;  Ct]) — — — , 

(sv)kn(Sv) 

u4>  =  U{0  — F<t)(sv){6o',  O',  <f>;  On)- 


R 


(2.20) 


(2.21) 


(2.22) 

(2.23) 

(2.24) 


where 


'R(SV)  = 


<5C33  —  SC  13 
P 


sin  0o  cos  0Q  cos2  0  — 


SC'i  i  —  SC  i3 


sin  0O  cos  6 o  sin2  0 


T  ~-6--  sin  0o  cos  0o  sin2  0  cos2  <b  -  cos  2 0o  sin  20  sin  <f> 

H  [i 


55 


*8(SV)  = 


8C33  —  8C\3 


9 


sin  0O  cos  0o  sin  9  cos  9  -f 


8CU  -  8C\Z 


9 


66  sin  0O  cos  $0  sin  29  cos2  <j>  H — cos  20o  cos  20  sin  <fr 
9  9 


sin  0O  cos  0o  sin  0  cos  0 

(2.26) 


F$(sv)  =  — sin  0O  cos  0q  sin  0  sin  2<£ -t — cos  290  cos  0  cos  <j> 


9 


9 


(2.27) 


Figures  7-9  display  three-dimensional  amplitude  patterns  resulting  from  the  interaction  of 
the  obstacle  with  the  impinging  SV  wave:  Fr(sv){90  —  0)  in  Figure  7;  Fq(sv){9q  =  0)  in 
Figure  8;  /\/,(svr)(0 0  =  0)  in  Figure  9.  Numerical  values  of  the  participating  parameters  are 
shown  in  Table  1. 


(c)  Incident  SH-wave  (fi  =  y/9/p) 


We  take 


7(u) 


—  iw[f- ^(ysinSo  +  ^cosPo)]  _  (SH)- 

—  t'o  cre  1  J  —  UQ  tj 


Since  the  non- vanishing  derivatives  are 


‘M'  -.(SI/),  ■  n  di 40) 

‘m3  SHI  0O  ;  — — 


(s'//) ,  n 
=  -?»o  0.3  cos  0O 


(2.28) 


(2.29) 
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only  and  in  (2.7)-(2.8),  survive.  We  then  find  that  in  this  case 


m  =  ?-FHsmV«  «;  c.v)^  (2.30) 

47 TO  it 

=  _t/(s//)^M^)^(sw)(^o;  ^  Ctj)£_^.  (2.31) 

=  (^o;  0;  (2.32) 

<i7T  IX 


with 


Fr{sh)  = 

<sc44 

cos  0O  sin  20  cos  v>  H - sin  0O  sin2  0  sin  2<£ 

(2.33) 

// 

Fg(SH)  = 

£C44 

0C66 

cos  0o  cos  20  cos  <f>  H - -  sin  0O  cos  0  sin  0  sin  2<£ 

(2.34) 

M 

/* 

F<t>(SH)  = 

6C44 

cos  0O  cos  0  sin  <f> - —  sin  0O  sin  0  cos  2  4> 

(2.35) 

Figures  10-12  show  the  radiation  pattern  resulting  from  the  scattering  of  SH  waves  by  the 
obstacle:  FR(Sh){6 o  =  90°)  in  Figure  10;  F9(SM)(0o  =  90°)  in  Figure  11;  F^SH)(0o  =  90°)  in 
Figure  12. 


3.  EVALUATION  OF  THE  ANISOTROPY  COEFFICIENTS 


The  scattered  field  is  governed  by  the  nine  elements  of  the  scattering  matrix 

Frp  Frsv  Frsh 

Fep  Fesv  Fbsh 
F<t,p  F^sv  Fj,sh 

These  elements  are  functions  of  five  independent  material  constants, 

6\  —  SC13  , 

82  —  8C33  ~  8C\$  , 

63  =  oCu  —  8C\z  »  (3-2) 

84  =  8C 44  , 

8 5  =  8Cee  ■ 

We  shall  prescribe  six  hypothetical  experiments  through  which  6Cn,  8C33,  8C13,  8C44,  8Cec 
arid  0o  can  be  determined. 
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Experiment  1:  Determination  of  the  Symmetry  Axis 

We  radiate  the  inclusion  with  a  plane  P-wave  and  observe  the  field  along  the  direction  of 
the  wave  normal  p  only,  i.e.:  at  the  point  where  9  =  ±60.  It  is  clear  from  Eqs.  (2.1 7)— (2. 19) 
that  as  we  vary  p  on  a  unit  sphere,  there  will  be  only  one  direction  in  space  along  which 
Fgp  —  F^p  —  0,  Fpp  ^  0.  This  direction  coincides  with  the  symmetry  axis,  relative  to  which 
9q  and  6  are  measured  and  {5V,  SH}  are  defined. 


Experiment  2:  Determination  of  SC44 


Radiate  the  inclusion  with  a  plane  5E  wave  along  Qq  =  0  and  observe  the  field  at  6  —  60  —  0. 
Then  via  (2.25)-(2.27) 


r  A  2 6C44  .  ,  n  2 8C44 

r  RSV  =  V  1  r$SV  — - sin  <p  \  r  <t>sv  =  - COS  0, 

p  p 


(3.3) 


will  yield  ( SC44 )  for  a  given  4>. 


Experiment  3:  Determination  of  6C&e 

Radiate  the  inclusion  with  a  plane  SH  wave  along  60  =  90°  and  observe  at  0  =  90°.  Then 
via  ( 2.33)— (2.35 )  we  find 

F rsi!  =  ——sin  2^  ;  F$su  =  0  ,  F^sn  = - —  cos  2d-  (3.4) 

p  p 
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Experiment  4:  Determination  of  6C13 


Radiate  the  inclusion  with  a  plane  P-wave  at  0O  =  0  and  observe  the  field  at  9  =  90°.  Then, 
via  (2.1 7)— (2. 19)  we  have 

FRp=irk  ’  fsp=0  •  F<p=0' 

Experiment  5:  Determination  of  6C33 

Radiate  the  inclusion  with  a  plane  P-wave  at  60  =  0  and  observe  the  field  at  6  =  0.  Then, 
via  (2.17)— (2.19): 

Frp  =  xtI;  ;  F,p"°  '• 

Experiment  6:  Determination  of  6Cu 

Radiate  the  inclusion  with  a  plane  P-wave  at  0Q  =  90°  and  observe  at  0  =  90°</>  =  90°.  Then, 
via  (2. 1 7)— (2. 19): 

Frp  -  ~2~  ;  Ftp  =  0  ;  F^r  =  0. 


60 


DISCUSSION 


The  theoretical  results  discussed  above  show  the  effect  on  elastic  wave  propagation  of 
Rayleigh  scattering  due  to  azimuthally  isotropic  obstacles  in  an  otherwise  homogeneous 
medium.  This  effect  may  be  important  for  the  study  of  propagation  of  elastic  waves  in  vari¬ 
ous  materials.  For  example,  the  differences  between  the  effects  of  azimuthally  isotropic  and 
isotropic  scatterers  can  be  considered.  The  radiation  pattern  corresponding  to  an  isotropic 
obstacle  will  be  the  same  as  that  presented  in  Figs.  2  and  4  for  vertically  incident  P  and 
SV  signals  on  the  azimuthally  isot.opic  scatterer.  As  the  angle  of  incidence  0 0  changes,  the 
isotropic  radiation  patterns  for  scattered  P  and  total  S-waves  will  simply  rotate  along  with 
the  incident  wave  number  vector,  although  there  will  be  some  interchange  between  SH  and 
SV  polarizations.  Comparison  of  Figs.  2  and  3  shows  that  the  radiation  patterns  for  the 
azimuthally  isotropic  scatterer  will  actually  change  shape,  depending  on  the  relative  mag¬ 
nitudes  of  the  different  perturbations  to  the  elastic  constants  bC\j.  The  most  important 
difference,  however,  is  that  since  the  scattering  amplitudes  for  horizontal  incidence  depend 
on  different  elastic  constants  than  for  vertical  incidence,  the  magnitude  of  backscattering 
may  be  significantly  different  for  these  two  cases.  On  the  other  hand,  for  the  isotropic  case, 
the  amplitude  will  remain  the  same. 

A  second  area  of  importance  for  the  scattering  theory  is  in  reference  to  theories  for 
inversion  of  acoustic  and  seismic  reflection  data.  An  analysis  of  radiation  patterns  for  point 
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difTractors  by  Tarantola  (1986)  shows  that  for  an  isotropic  medium,  an  earth  model  is  best 
parameterized  in  terms  of  S  and  P  impedances  and  density  in  order  to  perform  a  general 
inversion  for  medium  properties  using  nearly  vertically  propagating  waves.  Impedance  for 
an  anisotropic  medium  is  a  somewhat  ill  defined  concept  due  to  the  directional  variation  of 
velocity,  but  if  the  impedance  is  defined  using  vertical  velocities,  the  conclusions  of  Tarantola 
(1986)  can  still  be  applied  to  the  azimuthally  isotropic  medium.  The  constant  perturbations 
SC33  and  <5C44  would  ,ake  the  roles  of  6\  and  hp,  respectively,  in  this  case,  since  these  are 
the  constants  affecting  vertical  propagation. 

However,  a  determination  of  all  5  elastic  parameters  is  desirable.  While  the  list  of  hy¬ 
pothetical  experiments  in  Section  3  shows  that  in  principle  all  5  constants  are  attainable 
through  observations  of  scattering,  it  is  also  clear  that  a  wide  range  of  incidence  angles 
relative  to  the  axis  of  symmetry  and  three  types  of  incident  waves  are  required  for  the  inves¬ 
tigation.  Therefore,  experiment  within  the  earth  must  consider  both  nearly  vertically  and 
nearly  horizontally  propagating  signals,  and  a  large  data  set  must  be  obtained. 

Although  the  application  of  these  scattering  concepts  to  inversion  work  will  be  difficult, 
they  will  still  be  important  in  understanding  observations  and  forward  modeling  of  propa¬ 
gation  within  a  medium  containing  azimuthally  isotropic  inclusions,  such  as  shale  materials. 

1  he  variation  of  the  radiation  patterns  from  the  isotropic  case  may  help  to  explain  some 
unusual  observations,  such  as  anomalous  S-wave  signals  in  explosion  data  (e.g.,  Kisslinger 


1 1  al,  1961).  Considering  the  prevalence  of  anisotropic  media  in  many  areas  (Winterstein. 
1986),  inclusion  of  these  effects  is  warranted. 
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Table  I.  Parameter  values  used  to  calculate  the  radiation  patterns  presented  in 


Figures  2-12. 


Parameter 

Value 

Dimensions 

^0 

22.05  x  1010 

dyn  •  cm-2 

P  o 

35.97  x  1010 

dyn  •  cm'2 

P 

2.70 

g  •  cm-3 

UqLH2Sv 

1.0  x  106 

4  -2 

cm  •  sec 

1.0  x  109 

X 

sc33 

2.0  x  10° 

X 

SC\3 

0.5  x  109 

dyn  ■  cm-2 

6C\i 

1.0  x  109 

X 

2.0  x  109 

X 

I 

I 


Figure  2:  Radiatioa  pattern  for  scattered  P  waves  generated  by  a  P-wave  incident  along  the  v.m- 
metry  axis  {Oq  =  0°)  for  the  parameter  values  given  in  Table  1.  The  upper  portion  of  the  linn  re 
shows  a  vertical  cross-section  of  the  radiation  pattern  in  the  x  —  ;  plane,  while  the  lower  portion 
shows  a  perspective  of  the  radiation  in  three  dimensions.  For  this  incidence  angle,  the  radiat  ,ou 
has  rotational  symmetry  about  the  tr-axis. 
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Figure  3:  Radiation  pattern  for  '.-altered  SY- waves  generated  by  a  P-wave  incident  along  the 
symmetry  axis  (#o  —  05h  I  !.<■  medium  parameters  are  given  in  Table  1.  The  tipper  portion 
shows  a  vertical  cross-section  of  tie-  radiation  pattern,  and  the  lower  portion  shows  a  perspective 
of  the  radiation  in  three  dimensions. 
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Figure  5:  Radiation  pattern  for  scattered  v  waves  due  to  a  horizontally  incident  P-wave  ( 0t 
for  the  parameter  values  given  m  I  able  '  1  he  format  is  the  same  as  for  Fig.  4. 
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figure  10.  Radiation  pattern  for  scattered  P-w^ves  due  to  a  horizontally  incident  SH-wave 
JO  )  for  the  parameter  values  given  in  Table  J.  I  he  format  is  the  same  as  for  Fig.  4. 
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due  to  a  horizontally  incident  S 1 1  - 


Figure  11:  Radiation  na-’.-n  for  -  attend  SV- waves  du 

tgn  =  90°  i  for  the  paramet-r  '..ilu-s  m-n  in  l'ahle  1.  The  format  is  the  same  as  for  Fig.  » 


Figure  12:  Radiation  pattern  for  scattered  S II -waves  due  to  a  horizontally  incident  SH-wave 
(Qq  =  90°)  for  the  parameter  values  given  in  Table  1.  The  format  is  the  same  as  for  Fig-  4. 
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